LinG3D: visualizing the spatio-temporal dynamics of clonal evolution

Background Cancers are spatially heterogenous, thus their clonal evolution, especially following anti-cancer treatments, depends on where the mutated cells are located within the tumor tissue. For example, cells exposed to different concentrations of drugs, such as cells located near the vessels in contrast to those residing far from the vasculature, can undergo a different evolutionary path. However, classical representations of cell lineage trees do not account for this spatial component of emerging cancer clones. Here, we propose routines to trace spatial and temporal clonal evolution in computer simulations of the tumor evolution models. Results The LinG3D (Lineage Graphs in 3D) is an open-source collection of routines (in MATLAB, Python, and R) that enables spatio-temporal visualization of clonal evolution in a two-dimensional tumor slice from computer simulations of the tumor evolution models. These routines draw traces of tumor clones in both time and space, and may include a projection of a selected microenvironmental factor, such as the drug or oxygen distribution within the tumor, if such a microenvironmental factor is used in the tumor evolution model. The utility of LinG3D has been demonstrated through examples of simulated tumors with different number of clones and, additionally, in experimental colony growth assay. Conclusions This routine package extends the classical lineage trees, that show cellular clone relationships in time, by adding the space component to show the locations of cellular clones within the 2D tumor tissue patch from computer simulations of tumor evolution models. Supplementary Information The online version contains supplementary material available at 10.1186/s12859-024-05813-7.


Fig.S1. 3D lineage trees of individual clones drawn with MATLAB routines.
For each clone (denoted by a different color), the top row shows the 3D lineage tree with all cells belonging to that clone (LinG3DClone.mroutine), while the bottom row includes only the cells that survived to the end of the simulation (LinG3DAliveClone.m routine).A-A': initial clone #0; B-B': mutated clone #1; C-C': mutated clone #2; D-D': mutated clone #3; E-E': mutated clone #4; F-F': mutated clone #5. .

Supplemental Material S2.
Examples of 3D lineage trees in MATLAB, Python, and R from a simulation of a growing tumor with a large number of cellular clones (probability of mutation:   =0.05).
Since this simulation has generated a large number of clones (147), we focus here only on a small subset (12) of the largest clones and present their 3D lineage trees.Fig. S4-S6 show pairs of 3D lineage trees using our LinG3D routines in MATLAB (Fig. S4), Python (Fig. S5), and R (Fig. S6).Top rows show tress that include traces of all cells belonging to a given clone (Fig. S4-S6A-L); bottom rows show trees that include only traces of cells that survived to the end of the simulation (Fig. S4-S6A'-L').As in the previous case, each clone is denoted by a different color.These clones are smaller and much less dispersed than in the previous example since the cells have mutated more often and generated more clones.Again, the trees with alive cells are smaller than the trees with all cells, since many cells have died (due to random cell death or drug-induced cell death) or have left the computational domain.In some cases, there are no surviving cells, including the initial drug-sensitive clone (Fig. S4-S6A').The late clones, that is the clones emerging later during the tumor development, are smaller, and the majority of clones have not survived to the end of the simulation-more than 80 out of 147 (results not shown).

Examples of cellular clones and 3D lineage trees from cell culture
To illustrate the use of LinG3D routines for experimental data, we traced ten cellular clones in the series of bright field images of the 2D culture of U2OS sarcoma cells.The images were acquired automatically every 3 minutes for 96 hours, and eight annotated images from this experiment are shown in Fig. S7A-H.For each annotated image, we recorded coordinates of all cells arising from the 10 selected precursor cells (shown in Fig. S7A) and coordinated of two daughter cells if the traced cell has divided.The full 3D lineage trees for nine clones are shown in Fig. S7I-R.Their colors correspond to the color of the selected cells.

Mathematical equations from a model for S1 and S2
To illustrate how the 3D lineage tree algorithms can be applied, we developed a relatively simple agent-based model of the growing tumor exposed to a cytotoxic drug (MATLAB routine: tumorGrowth_example1.m and tumorGrowth_example2.m).This model is based on our Multi-Cell Lattice-Free, MultiCell-LF, framework [1][2][3] in which individual cells interact with one another through physical forces and with the surrounding microenvironment via chemical factors, such as drugs or oxygen.

Computational domain
The model is defined on a square tissue patch: [−100,100] × [−100,100]  $ , that contains five irregularly distributed stationary vessels  % ( = 1, … ,5) that are the source of a cytotoxic drug (, , ).The tumor cells  % ( = 1, … ,  ' ) can proliferate, absorb the drug, undergo random and drug-induced death, and can mutate which results in cell resistance to the drug.The initial condition consists of a single tumor cell located in the middle of the domain with no drug.The noflux boundary conditions are imposed on the drug domain.The cells are also removed from the system if they move outside of the domain boundaries.

Drug kinetics
Drug distribution within the tumor tissue depends on three factors: the amount of drug supplied by the vessels, the amount of drug taken up by the tumor cells, and the spatial localization of all blood vessels and cells.The drug (, , ) is supplied from each vessel  % with the influx rate ℐ ( , it diffuses through the tissue with a constant diffusion coefficient  ( , and is absorbed by each cell  % with the uptake rate  ( proportional to the drug concentration (, , ).The equation of drug kinetics is defined as follows: where interactions between the drug grid (, ) and the individual cells or vessels ( =  % or  % ) are specified by the indicator function with radius ℛ:

Cell division and mutation
Each cell is defined by its position  % and a constant radius ℛ ' .The cell can inspect its vicinity of radius ℛ ) and can count the number  % of neighboring cells.Cell age progresses with time ( * !*, = 1), and when the cell reaches maturity,  % ≥  % *%-, and is not overcrowded, that is the number of cell neighbors does not exceed the prescribed threshold,  % ≤  ./0, the cell will divide.Upon division of cell  % , two daughter cells  % " and  % # are created. % " will take position of the mother cell, i.e.  % " =  % and  % # will be located in the vicinity of the mother cell in the randomly selected direction, i.e.,  % # =  % + ½ ℛ ' (, ), where  ∈ [0,2].
The mutation status of each daughter cell is inherited from the mother cell, i.e., ℳ % " = ℳ % # = ℳ % .However, each daughter cell can acquire a new mutation status when it is born; the probability of new mutation is determined by  .1,(the only parameter we varied in the presented examples).Cell mutation will result in a shorter cell division age, i.e.,  % *%-= ½  ./,, where  ./, is the default cell division (maturation) age.It is assumed that the random mutations will start after the cell colony reaches a noticeable size.

Cell drug absorption, accumulation, and drug-induced death
Each cell absorbs drug (, , ) from its immediate vicinity of radius ℛ ( .The cell uptake rate is  ( .Drug accumulation  % by the cell  % is described by the following equation: • (, , )  ℛ # ((, ),  # ) Upon cell division, the levels of absorbed drug for both daughter cells are set to 0, i.e.,  % " =  % # = 0.When the drug accumulated by a non-mutated cell (ℳ % = 0) exceeds the prescribed threshold, i.e.  % >  ,23 , the cell will die with a probability   *%5 .Once the cell dies, it is removed from the system.

Cell random death
All cells, mutated and non-mutated, can undergo random death with a probability  *%5 , provided that they are of considerable age, that is their current age  % exceed twice the cell division age  % *%-.Once the cell dies, it is removed from the system.

Cell-cell and cell-vessel interactions
Individual cells  % and  6 can exert a cell-cell repulsive force  %,6 (with stiffness  ' and a resting length 2ℛ ' ) to avoid overlapping: If the cell  % overlaps with multiple (M) cells, the cumulative repulsive force is: To avoid overlapping between cells and vessels, a repulsive force  %,6 (with stiffness  -and a resting length (ℛ ' + ℛ -)) between the overlapping cell  % and the vessel,  6 is exerted on the cell: If cell  % overlaps with multiple (M) vessels, the cumulative repulsive force is

Cell Relocation
Cell dynamics is governed by Newton's second law of motion, where the applied forces arise as a result of cell-cell repulsive interactions  % 358 , repulsive interactions between cells and vessels  % 358 , and forces needed to overcome medium viscosity  % -%9 , where  % is cell mass: Assuming that springs are overdamped:  % * # : !*, # = 0, and that viscous force is proportional to cell velocity:  % -%9 = − *: !*, , the equation of cell relocation is given by: All model parameters are listed in Table 1.Output data used by the 3D lineage tree algorithms In order to draw the 3D lineage trees, the output data needs to be prepared in the following format: cell_history=[cell ID, clone ID, mother ID, birth iter, div/death iter] cellXY_iter=[x, y] cellID_iter=[cell ID] and saved as text files: cell_history.txt,cellXY_iter.txt,and cellID_iter.txt,respectively, where iter is the iteration number.Here, cell ID is a unique ID number for the cell, clone ID is a number unique to a given clone to which the cell belongs, mother ID is a unique ID number of the cell's mother cell, birth iter is an iteration number at which the cell was born, div/death iter is the iteration number at which the cell either divided into two daughter cells or died, x and y are the coordinates of the cells at a given iteration iter.The vector cell_history contains history of all cells from the whole simulation and should be save at the end of simulation.The vectors cellXY_iter and cellID_iter are created and saved for the iteration with number iter.These data sets are initiated by: cell_history=[1,0,0,0,0], which means that the first cell has unique index 1, belongs to clone number 0, cell's mother index is 0, the cell was born in iteration 0, and is still alive (has neither died nor divided); cellXY_0=[x0,y0] are the coordinates of the initial cell; cellID_0= [1] is a unique index of the initial cell (first cell).Subsequently, these data sets have to be updated in the following situations: 1. Cell division: (a) for the mother cell, find the row in cell_history that corresponds to the unique mother ID and update the last element to indicate the current iteration number, i.e., iteration at which the cell has divided; (b) for each daughter cell, create a new row in cell_history and add new unique cell ID (keep track of the last unique ID used), copy clone ID from the mother cell, copy mother cell unique ID, add the current iteration number to indicate iteration at which the cell was born, keep last element equal to 0 to indicate that the cell is alive.Repeat the same for the second daughter cell (cell ID must be unique).2. Cell mutation: find the row in cell_history that corresponds to the unique cell ID and update the second component with a new unique clone ID number (keep track of the last unique clone ID used).3. Cell random death: find the row in cell_history that corresponds to the unique cell ID and update the last component to indicate the iteration in which the cell died.4. Cell drug-induced death: find the row in cell_history that corresponds to the unique cell ID and update the last component to indicate the iteration in which the cell died. 5. Cell removal from the system if it moved outside the computational domain: find the row in cell_history that corresponds to the unique cell ID and update the last component to indicate iteration in which the cell was removed.6. Saving cell data: for those iterations iter, for which data is being saved in the external text files, save coordinates of all cells currently in the system (file cellXY_iter.txt)and unique cell IDs (file cellID_iter.txt) in the very same order in which the coordinates are saved.7. Saving drug data: for the last iteration, save the values of drug concentration in all grid points (in (external text file drug.txt).